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1  Introduction 


Observations  have  shown  that  convecting  storms  create  significant  internal  wave  ac¬ 
tivity;  see  for  example  the  recent  satelite  observations  of  Dewan,  et  al  [1].  Direct 
numerical  simulations  [2,  3,  4,  5]  also  have  identified  convective  storms  as  a  source 
of  gravity  waves.  There  is  little  doubt  that  convection  creates  wave  energy  in  the 
atmosphere.  Internal  waves  in  the  atmosphere  may  be  created  by  other  phenomena, 
as  well  as  convection.  The  other  phenomena  most  often  blamed  for  wavemaking  are 
jet  streams,  fronts,  and  flow  over  mountains.  The  relative  importance  of  the  different 
wave  sources  is  a  topic  of  current  debate. 

The  mechanism  of  gravity  wave  generation  by  convecting  storms  is  not  clear.  Clark 
[2]  claims  there  are  two  mechanisms.  The  first  mechanism  concerns  the  zonal  flow 
(average  horizontal  flow)  that  usually  exists  in  the  atmosphere.  Convection  acts  as 
an  obstruction  to  this  flow,  deflecting  the  airflow  over  and  around  the  thunderstorm, 
analogous  to  flow  over  a  mountain.  This  obstruction  creates  a  pattern  of  steady 
waves,  analogous  to  mountain  waves. 

The  second  mechanism  concerns  an  oscillation  of  the  flow  inside  a  thunderstorm, 
which  can  act  as  a  mechanical  wavemaker.  Fovell,  et  al,  [4]  studied  the  second  mecha¬ 
nism  using  an  artificial  mechanical  oscillator  in  a  realistic  atmosphere.  The  results  of 
Fovell  are  two-dimensional  direct  numerical  simulations.  Fovell  found  that  the  simple 
mechanical  oscillator  created  a  fixed  V-shaped  pattern  of  waves.  Fovell  showed  that 
the  wave  pattern  is  largely  influenced  by  the  zonal  shear  flow,  and  by  the  orienta¬ 
tion  of  the  mechanical  oscillator.  They  argue  that  the  stratospheric  waves  apparent 
in  the  simulations  are  the  result  of  an  interaction  between  the  convection  and  the 
tropopause.  Realistic  simulations  of  convecting  storms  have  also  been  performed  by 
Lane,  et  al  [3],  and  Holton  and  Alexander  [5].  The  V-shaped  wave  pattern  is  clearly 
evident  in  both  sets  of  results,  similar  to  the  waves  produced  by  Fovell. 

Numerical  simulations  of  a  simple  model  of  convection  are  considered  below.  The 
simulations  are  restricted  to  two  spatial  dimensions.  The  convection  is  modeled  as  a 
vortex  pair,  scaled  to  match  the  typical  size  and  strength  of  a  real  convecting  storm. 
However,  the  tremendous  complexity  of  an  actual  storm  is  not  included.  Instead,  by 


including  only  the  simple  model  of  the  large  scale  flow,  and  the  resulting  gravity  wave 
activity,  the  relationship  between  convection  and  waves  is  more  easily  understood. 

Previous  numerical  simulations  by  Hill  [6]  also  considered  a  vortex  pair  in  a  strat¬ 
ified  fluid.  Waves  are  not  mentioned  in  Hill’s  work,  and  the  results  focus  on  the 
movement  of  the  center  of  the  vortex.  The  size  and  strength  of  the  vortices  of  Hill 
were  chosen  to  represent  the  motion  of  trailing  vortices  behind  an  aircraft,  rather 
than  the  large  scale  motion  of  storms.  The  stratification  is  weak  when  considering  an 
aircraft’s  trailing  vortex.  The  scale  of  a  convecting  storm  is  large,  and  the  stratifica¬ 
tion  is  relatively  strong,  which  accounts  for  the  difference  between  Hill’s  results  and 
the  results  given  below. 

The  equations  treated  here  are  the  anelastic  equations.  The  anelastic  equations 
model  the  dynamics  of  an  adiabatic  atmosphere  approximately,  and  do  not  allow 
sound  wave  propagation.  The  lack  of  sound  waves  does  not  seriously  affect  the  dy¬ 
namics  of  the  internal  waves,  but  does  allow  a  significantly  larger  time  step  when 
compared  to  a  numerical  model  of  the  fully  compressible  equations.  The  anelastic 
equations  were  first  suggested  by  Batchelor  [7],  and  then  further  developed  by  several 
authors  [8], [9], [10].  The  original  anelastic  equations  have  a  flaw;  the  linear  solution 
does  not  match  the  linear  solution  of  the  fully  compressible  equations.  It  was  later 
shown  that  the  pressure  term  was  incorrect.  Bacmeister  and  Schoeberl  [11]  provide 
a  concise  derivation  of  the  corrected  anelastic  equations. 

The  numerical  method  is  spectral  in  space  with  periodic  sidewall  boundaries,  and 
rigid  boundaries  on  the  top  and  bottom.  A  wave  damping  layer  is  also  included  near 
the  top  boundary.  The  temporal  algorithm  is  semi-implicit,  and  uses  a  formulation 
where  the  pressure  and  the  horizontal  velocity  are  eliminated  from  the  linear  terms 
in  the  equations.  Details  are  provided  by  McHugh  [12]. 

Three  base  states  are  considered;  a  constant  density  case,  a  constant  Brunt- Vaisala 
case,  and  a  case  with  two  layers,  each  with  its  own  constant  value  of  the  Brunt- Vaisala 
frequency.  This  last  case  models  the  atmosphere  in  the  region  of  the  tropopause. 


2  Governing  equations 


Although  the  results  given  later  are  two  dimensional,  the  equations  and  numerical 
methods  are  capable  of  three  dimensions,  and  are  developed  here  in  three  dimensions. 
The  reduction  to  two  dimensions  is  discussed  where  appropriate. 

The  anelastic  equations  for  a  compressible  atmosphere  are 
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u,  v,  w  are  the  velocity  components,  x,  y,  z  are  the  components  of  position,  0  is  the 

potential  temperature,  p  and  0  are  the  basic  state  density  and  potential  temperature, 

Cp  is  the  specific  heat  at  constant  pressure,  R  is  the  gas  constant,  p  is  the  pressure, 

and  po  is  a  constant.  Equations  (1-3)  are  the  momentum  equations,  (4)  is  the  energy 

equation,  and  (5)  is  the  continuity  equation.  They  have  been  rescaled  using  a  velocity 

scale,  U,  to  be  defined  later,  and  a  length  scale,  L,  which  is  the  horizontal  length 

(parallel  to  the  x  direction)  of  the  domain.  Note  that  Lx,  Ly,  and  Lz  are  length  scale 


ratios  which  account  for  the  difference  between  the  physical  domain  lengths  and  the 
computational  domain  lengths.  The  density  and  potential  temperature  are  rescaled 
using  their  respective  values  at  the  bottom  of  the  domain  (po,  90).  The  Reynolds, 
Prandtl,  and  Froude  numbers  are 


& 

n> 

II 

J  ^ 

(10) 

^  fse 
II 

(11) 

U2 

F2  =  — . 
r  gL 

The  basic  state  is  governed  by  the  perfect  gas  law, 

(12) 

p  =  pRT , 

(13) 

and  the  equation  of  static  equilibrium, 

dp 

Yz  =  ~~p9- 

(14) 

Generally  the  base  state  temperature  profile,  T,  is  chosen,  and  the  remaining  base 
state  variables  are  determined  using  (7),  (13),  and  (14). 

One  important  base  state  parameter  is  the  Brunt- Vaisala  frequency,  defined  as 


9 dQ_ 

9  dz' 

Rescaling  this  frequency  as  before  gives  the  Brunt- Vaisala  parameter: 

2=L£l^  =  _^ia0 

U29dz  F2  9  dz  ’ 

where  z  is  now  dimensionless.  Another  base  state  parameter  is  the  scale  height, 
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3  The  Numerical  Scheme 


The  numerical  method  is  that  of  McHugh  [12].  It  is  similar  to  the  normal  velocity- 
normal  vorticity  method  of  Kim,  Moin,  and  Moser  [13],  who  studied  the  incompress¬ 
ible  Navier-Stokes  equations.  The  governing  equations  are  reduced  such  that  pressure 
and  horizontal  velocity  components  are  eliminated  from  the  linear  terms,  resulting  in 
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where  Ai  is  the  sum  of  the  nonlinear  terms  for  the  ittl  momentum  equation, 
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The  variable  <j)  is  retained  in  the  calculations,  and  w,  <j>,  and  9  are 

determined  si- 

multaneously  using  (21),  (22),  and  the  energy  equation.  The  energy  equation  is  now 

written  as 
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where  B  is  the  sum  of  the  nonlinear  terms  for  the  energy  equation, 
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The  boundary  conditions  are 

w  =  wz  =  9  =  0 

(25) 
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on  the  rigid  boundaries.  Note  that  the  sidewalls  are  treated  with  periodic  boundary 
conditions,  and  no  further  enforcement  is  allowed. 

After  m,  <j>,  and  9  are  found,  the  remaining  task  is  to  determine  the  horizontal 
velocity,  u  and  v.  There  is  a  particularly  straightforward  method  of  finding  the 
horizontal  velocity  in  two  dimensions.  Keeping  x  as  the  horizontal  direction  and  z  as 
the  vertical,  the  continuity  equation  is 


This  equation  is  easily  solved  for  u  once  w  is  available. 

In  three  dimensions,  the  vorticity  equation  is  necessary,  along  with  the  continuity 
equation,  to  determine  both  components  of  velocity.  The  vertical  vorticity  equation 
for  viscous  flow  is 
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where  77  is  the  vertical  component  of  vorticity.  The  boundary  condition 


on  the  rigid 


boundaries  is 


77  =  0. 


(28) 


Equation  (27)  is  solved  for  the  vorticity,  subject  to  (28),  and  then  the  definition  of 
vorticity  provides  one  equation  for  horizontal  velocity: 


Uy  -Vx  =  7 7.  (29) 

This  result,  along  with  the  continuity  equation, 

(30) 

are  sufficient  to  determine  u  and  v. 

The  spatial  discretization  is  a  spectral  method.  The  horizontal  directions  are  ex¬ 
panded  in  Fourier  series.  The  vertical  direction  uses  a  spectral  element  method,  using 
Lagrange  interpolants  within  each  subdomain,  collocated  on  the  Chebyshev-Gauss- 
Lobatto  points.  A  damping  layer  is  used  at  the  top  of  the  domain  as  a  non-reflecting 
lid.  Note  that  the  damping  layer  is  ineffective  for  gravity  waves  with  wavelengths 
greater  than  the  damping  layer  thickness. 
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4  Initial  Conditions:  A  vortex  pair 


The  initial  conditions  are  a  pair  of  point  vortices.  The  velocity  field  for  a  point  vortex 
is  taken  from  the  well-known  similarity  solution  for  a  diffusing  line  vortex,  as  discussed 
in  Sherman  [14].  This  similarity  solution  is  an  exact  solution  to  the  Navier-Stokes 
equations  for  constant  density,  viscous  flow.  In  cylindrical  coordinates,  the  velocity 
field  for  a  single  vortex  located  at  the  origin  is 
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where  v  is  now  the  azimuthal  velocity,  r  is  the  radial  distance  from  the  center  of  the 
vortex,  T  is  the  circulation  around  the  vortex,  and  v  is  the  kinematic  viscosity.  The 
radial  component  of  velocity  is  zero.  This  solution  is  rescaled  using  the  same  length 
and  velocity  scale  as  before.  The  resulting  azimuthal  velocity  is 


v  = 


27rr 


1  —  t 


(32) 


where  v,  r,  and  t  are  now  dimensionless,  and  the  dimensionless  vortex  strength  is 
defined  as 


G  = 


T 

LU' 


(33) 


Note  that  the  solution  given  by  (32)  depends  on  time  as  well  as  space.  Only  one 
value  of  time  was  used  to  start  the  simulations.  The  vortex  was  merely  used  to  start 
a  flow  that  has  the  general  characteristics  of  the  vortex  pair  being  studied. :  : 

An  obvious  alternative  to  the  viscous  vortex  is  the  well-known  irrotational  inviscid 
vortex.  Similations  with  the  irrotational  vortex  led  to  numerical  instabilities,  clearly 
originating  at  the  core,  where  the  analytic  solution  is  unbounded.  Attempts  to  smooth 
the  initial  conditions  at  the  core,  and  then  advance  the  solution  with  the  chosen 
spectral  method,  were  tedious  and  only  marginally  successful.  Hence  the  irrotational 
vortex  was  not  pursued. 


5  Results  and  Discussion 

The  flow  is  initiated  by  releasing  a  pair  of  counterrotating  vortices  in  an  otherwise 
motionless  flow.  The  vortex  pair  is  chosen  to  rotate  with  a  downward  velocity  in  the 


Table  1:  Dimensionless  parameter  values 


Reynolds  number 

2000 

Prandtl  number 

1 

Froude  number 

1 

Vortex  strength 

1 

Lx/  L, 

1 

region  between  the  vortices.  This  choice  is  based  on  the  initiation  of  a  storm  event 
.in  the  atmosphere,  which  is  well-known  to  result  in  a  downward  wind  in  the  storm 
center.  The  vortices  are  initially  located  at  one  quarter  of  the  domain  height  from 
the  bottom  (z=-0.5),  and  a  dimensionless  distance  of  1.0  from  the  horizontal  center. 

There  is  no  zonal  flow  in  the  present  simulations,  which  means  that  a  velocity  scale 
must  be  selected  from  other  parameters.  The  velocity  scale  is  chosen  to  be  U  =  \/gL, 
where  g  is  the  gravitational  constant,  and  L  is  horizontal  domain  length.  The  Froude 
number  becomes  unity,  the  Prandtl  number  and  the  scale  height  are  unchanged,  and 
the  other  parameters  are  now 


vW 

’ 

If:  ;  (34) 
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vV'3’ 
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Before  discussing  the  results  for  the  stratified  case,  it  is  useful  to  consider  a  flow 
without  stratification.  The  base  state  density,  potential  temperature,  and  pressure- 
are  chosen  to  be  constant,  and  the  flow  is  started  with  the  same  initial  conditions 
as  the  stratified  case.  Since  the  boundary  conditions  on  temperature  are  homoge¬ 
neous,  and  the  initial  temperature  field  is  zero,  the  potential  temperature  remains 
zero  throughout  the  simulation  for  this  case.  The  total  effect  is  that  the  equations 
are  reduced  to  the  incompressible  Navier-Stokes  equations.  Table  1  lists  the  chosen 
dimensionless  parameter  values.  Table  2  lists  the  physical  parameter  values  on  which 
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Figure  5:  Contour  plot  of  vorticity  for  the  constant  density  case  with  t=1.0 

H0.20 

0.10 

0.00 

-0.10 

-0.20 


Figure  6:  Contour  plot  of  vorticity  for  the  constant  density  case  with  t=32.0 

■  0.20 

0.10 

0.00 

-0.10 

-0.20 


14 


Figure  7:  Base  state  density  profiles 


atively  short  time  span  of  t=32  dimensionless  units.  The  vortices  move  together 
slightly,  reminiscent  of  the  Bernoulli  effect  between  two  spinning  cylinders.  The  vor¬ 
tices  also  decay  in  strength,  and  diffuse  outward  from  the  center  of  the  vortex.  The 
close  proximity  of  the  two  vortices  causes  them  to  become  elongated  in  the  vertical, 
and  to  induce  an  upwards  motion  of  the  vortex  pair.  This  vertically  induced  motion 
of  a  vortex  pair  is  well-known,  and  is  discussed  in  Lamb  [15]. 

Now  consider  stratified  flow.  The  base  state  is  chosen  to  have  a  contant  Brunt- 
Vaisala  frequency,  achieved  with  a  linear  variation  in  temperature.  Figure  7  shows 
the  base  state  density  profile.  The  dimensionless  parameter  values  are  the  same  as 
those  for  the  constant  density  case,  as  shown  in  Table  2,  with  the  addition  of  the 
Brunt- Vaisala  parameter  value  of  0.32. 

The  results  of  this  stratified  case  are  shown  in  figures  8-19,  using  the  same  com¬ 
bination  of  velocity  vector  plots,  surface  plots  of  the  velocity  magnitude,  and  colored 
contour  plots  of  the  vorticity. 

The  dynamics  of  the  vortices  in  stratified  flow  are  dramatically  different  than  the 
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Figure  8:  Velocity  vector  plot  for  the  constant  N  case  with  t=1.0 


Figure  10:  Velocity  vector  plot  for  the  constant  N  case  with  t=25.0 
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Figure  11:  Velocity  vector  plot  for  the  constant  N  case  with  t=32.0 
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Figure  15:  Surface  plot  of  velocity  magnitude  for  the  constant  N  case  with  t=32.0 
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Figure  16:  Contour  plot  of  vorticity  for  the  constant  N  case  with  t=1.0 
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Figure  18:  Contour  plot  of  vorticity  for  the  constant  N  case  with  t=25.0 


Figure  19:  Contour  plot  of  vorticity  for  the  constant  N  case  with  t=32.0 


constant  density  case.  As  soon  as  the  vortices  are  released  in  the  stratified  flow,  they 
very  rapidly  lose  virtually  all  of  the  kinetic  energy,  as  shown  in  figures  12-15.  The 
surface  plot  in  figure  12  clearly  shows  the  two  initial  vortices,  although  they  have 
decayed  somewhat  when  compared  to  the  surface  plot  for  the  constant  density  case 
at  the  same  time,  shown  in  figure  1.  However,  the  coherency  of  the  vortices  disappears 
very  quickly,  and  is  replaced  with  an  internal  wave  pattern,  as  shown  in  figures  13-15. 
This  conversion  of  energy  in  the  vorticies  into  wave  energy  occurs  in  approximately 
one  Brunt- Vaisala  period, 


(37) 


defined  as  the  time  for  one  oscillation  at  the  Brunt- Vaisala  frequency,  and  has  a 
dimensionless  value  of  approximately  9.8  for  the  chosen  Brunt- Vaisala  parameter  value 
of  0.32.  The  wave  motion  continues,  and  the  initial  vortex  pair  is  never  reformed. 
Interaction  with  the  damping  layer  at  the  top,  and  viscous  dissipation,  finally  cause 
the  wave  action  to  decay  to  zero. 

It  should  be  noted  that  the  motion,  initially  symmetric  about  the  vertical  center 
of  the  domain,  retains  this  symmetry  during  the  entire  simulation.  Clearly  the  waves 
propagate  rapidly  to  the  sidewalls,  and  interact  with  the  motion  on  the  Other  side. 
Since  the  sidewalls  are  periodic,  one  may  think  of  an  identical  motion  on  each  side  of 
the  domain,  evolving  simultaneously.  The  results  are  therefore  also  symmetric  about 
the  vertical  sidewall,  and  when  wave  energy  impinges  the  sidewall,  it  is  perfectly 
reflected  due  to  this  symmetry.  This  effect  does  not  model  the  real  atmosphere, 
although  the  solutions  are  accurate  and  satisfy  all  equations  and  boundary  conditions. 
The  results  presented  here  focus  on  the  motion  before  this  reflection  occurs. 

Figures  5-6  show  the  time  evolution  of  the  vorticity  for  the  constant  density  case, 
The  initial  vortices  retain  their  identity,  but  dissipate,  elongate,  and  move  together 
and  upwards.  Note  that  the  vorticity  field  is  antisymmetric,  and  the  dark  and  bright 
spots  are  due  to  the  counterrotating  nature  of  the  vortex  pair. 

Figure  16  shows  the  vorticity  contours  for  the  stratified  case  after  a  mere  ten 
time  steps,  and  clearly  the  vortices  are  still  coherent,  and  nearly  identical  to  figure 
5.  Figure  17  shows  two  streaks  of  vorticity  after  an  elapsed  dimensionless  time  of 


18,  corresponding  to  the  time  for  several  oscillations  at  the  Brunt- Vaisala  frequency. 
Also  evident  in  figure  17  are  several  other  dark  and  bright  spots  of  vorticity,  in  the 
vicinity  of  the  streaks.  The  initial  vortices  do  not  evolve  into  these  dark  and  light 
streaks,  as  one  might  assume.  Instead,  the  initial  vortices  are  now  the  small  patches 
of  vorticity  visible  on  the  outside  of  the  larger  streaks.  The  streaks  of  vorticity  have 
been  generated  by  the  appearance  of  an  internal  wave  motion,  beginning  to  dominate 
the  flow  pattern.  If  the  vorticity  patches  from  the  initial  conditions  are  tracked,  they 
are  found  to  disappear  within  a  time  equivalent  to  three  oscillations  at  the  Brunt- 
Vaisala  frequency.  The  motion  of  the  vortex  center  is  nearly  horizontal,  and  moving 
away  from  the  domain  center. 

Figures  18  and  19  are  later  time  slices  of  the  vorticity,  showing  regular  patterns 
of  positive  and  negative  vorticity.  These  patterns  indicate  wave  motion,  rather  than 
any  historical  result  of  the  initial  vortices.  It  is  interesting  to  note  in  figure  19  that 
two  intense  patches  of  vorticity  appear  near  the  upper  corners  of  the  domain,  ap¬ 
proximately  equal  in  strength  to  the  original  vortex.  These  patches  appear  as  the 
wave  motion  reaches  its  greatest  amplitude,  as  can  be  seen  in  figures  12-15,  which  are 
surface  plots  of  the  magnitude  of  the  velocity  at  times  corresponding  to  figures  16-19. 
This  wave  is  not  breaking,  meaning  the  vertical  gradient  of  potential  temperature  is 
a  stable  gradient,  and  the  flow  at  this  Reynolds  number  is  not  unstable.  Yet  one 
wonders  if  this  motion  at  other  parameter  values  could  result  in  a  Helmholtz  shear 
type  instability  near  the  wave  crest,  and  a  resulting  patch  of  turbulence. 

The  third  base  state  is  a  two  layer  case,  with  constant  Brunt- Vaisala  frequency 
in  each  layer.  The  interface  between  layers  is  located  at  the  vertical  center  of  the 
domain.  Figure  7  shows  the  base  state  density  profile  for  this  case,  as  well  as  the 
previous  case.  The  parameter  values  are  identical  to  the  second  case,  except  now  the 
upper  half  of  the  domain  has  a  Brunt- Vaisala  parameter  value  of  0.62,  corresponding 
to  a  Brunt- Vaisala  frquency  of  0.02  1/s,  typical  of  the  stratosphere. 

The  results  of  the  two-layer  case  are  shown  in  figures  20  through  31,  using  the  same 
time  values  and  graphics  as  the  constant  N  case.  The  dynamics  of  the  atmosphere 
for  this  two  layer  case  is  very  similar  to  the  single  layer  of  constant  N.  The  vortices 
lose  coherency  very  quickly  and  form  an  internal  wave  pattern.  The  amplitude  of 


the  waves  in  the  upper  layer  for  this  two  layer  case  are  generally  greater  than  the 
amplitude  in  the  one  layer  case.  This  increase  in  amplitude  is  due  primarily  to  the 
lower  density  at  the  higher  elevations,  rather  than  any  dynamic  effect  at  the  interface 
of  the  two  layers.  Note  that  the  waves  in  the  upper  layer  reach  a  maximum  amplitude 
more  quickly  compared  with  the  single  layer  case. 

Perhaps  the  most  significant  difference  is  evident  in  the  contour  plots,  figures  28- 
31.  Figure  28  shows  the  vorticity  after  a  mere  10  time  steps,  and  should  be  compared 
with  figure  16,  the  vorticity  contours  for  the  single  layer  case  at  the  same  time. 
The  vorticity  with  two  layers  at  this  early  stage  is  nearly  identical  to  the  single  layer, 
except  at  an  elevation  just  above  the  interface  between  the  two  layers.  A  narrow  strip 
with  high  levels  of  vorticity  appears,  as  evident  in  figure  28.  The  dynamics  in  figure 
28  are  not  yet  dominated  by  internal  waves,  but  the  sudden  change  in  Brunt- Vaisala 
frequency  is  clearly  generating  a  layer  of  motion  with  enhanced  vorticity. 

Figure  29  shows  that  this  layer  of  enhanced  vorticity  has  evolved  and  generated 
two  spots  of  strong  vorticity,  in  addition  to  the  spots  and  streaks  of  vorticity  that 
existed  with  a  single  layer  shown  in  figure  17.  Except-  for  this  new  layer  of  vorticity 
just  above  the  interface  of  the  two  layers,  the  vorticity  pattern  matches  the  single 
layer  case  rather  closely. 

Figures  30  and  31  show  the  vorticity  contours  for  the  time  values  corresponding  to 
the  single  layer  case,  shown  in  figures  18  and  19.  Below  the  interface,  the  dynamics  of 
the  two  layer  case  appear  nearly  identical  to  the  single  layer  case.  Above  the  interface, 
where  the  Brunt- Vaisala  parameter  is  much  different,  the  wave  motion  evolves  on  a 
different  time  scale.  But  the  basic  motion  is  similar.  This  is  evident  by  the  strong 
spots  of  vorticity,  appearing  at  a  time  of  32  for  the  single  layer  case,  but  near  a  time 
of  25  for  the  two  layer  case.  In  fact  the  basic  patterns  of  vorticity  above  the  interface 
between  figure  19  and  figure  30  compare  rather  well,  considering  the  time  difference 
between  the  two  figures. 

The  dispersion  of  energy  from  the  original  vortices  is  difficult  to  quantify,  since  the 
vortices  become  so  ill-defined  in  a  very  short  time.  However,  the  flux  of  energy  across 
a  horizontal  plane  is  a  quantity  indicative  of  the  flux  from  the  original  vortex,  and 
provides  information  concerning  the  later  dynamics  of  the  waves.  Figure  39  shows  the 


Figure  20:  Velocity  vector  plot  for  the  two  layer  case  with  t=1.0 


Figure  24:  Surface  plot  of  velocity  magnitude  for  the  two  layer  case  with  t=1.0 


Figure  25:  Surface  plot  of  velocity  magnitude  for  the  two  layer  case  with  t=18.0 
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Figure  26:  Surface  plot  of  velocity  magnitude  for  the  two  layer  case  with  t=25.0 


28 


29 


3i 


temporal  evolution  of  the  vertical  energy  flux  for  all  three  base  states,  where  vertical 
energy  flux  is  defined  as 

Ew  =  w  ( u 2  +  w2),  (38) 

and  the  overline  represents  a  horizontal  average.  The  value  shown  in  figure  39  corre¬ 
sponds  to  the  energy  flux  across  the  vertical  center  of  the  domain. 

The  value  of  Ew  for  the  constant  density  case  starts  out  negative,  as  shown  in 
figure  39,  due  to  the  strong  downdraft  in  the  region  between  the  two  vortices.  As 
the  vortices  move  upward,  this  energy  flux  gradually  becomes  positive  and  then  stops 
increasing  as  the  position  of  the  center  of  the  vortices  approaches  the  vertical  center 
of  the  domain. 

The  value  of  Ew  for  the  stratified  cases  are  also  shown  in  figure  38.  Note  the 
dramatic  difference  between  the  constant  density  case  and  the  stratified  cases.  The 
single  layer  case  very  quickly  shows  large  positive  values  of  Ew,  associated  with  the 
rapid  transfer  of  vortical  energy  to  wave  energy.  Then  a  very  large  negative  energy 
flux  occurs.  This  negative  energy  flux  is  associated  with  the  complicated :  wave  mo¬ 
tion  in  the  region  between  where  the  original  vortices  were  located.  This  effect  occurs  _ 
before  the  elapsed  time  has  reached  even  one  Brunt-  Vaisala  period,  and  before  any. 
significant  energy  has  reached  the  boundaries.  The  remainder  of  the  data  shows  some 
oscillations,  centered  around  zero  energy  flux.  The  flow  beyond  a  time  of  approx¬ 
imately  32  is  likely  to  be  affected  by  reflected  waves  off  the  sidewalls.  Hence,  the 
later  wiggles  in  figure  38  may  not  be  indicative  of  motion  in  a  real  nonsymmetric 
atmosphere.  The  energy  flux  for  the  two  layer  case  is  very  similar  to  the  single  layer 
case.  The  rate  of  decay  of  the  energy  flux  as  the  initial  positive  value  changes  to  a 
negative  value  is  nearly  identical. 

Another  measure  of  energy  transfer  is  the  temporal  evolution  of  the  total  kinetic 
energy,  defined  as 

E  =  Jvp(u2  +  v2)dV,  (39) 

and  is  shown  in  figure  39,  ‘again  for  all  three  base  states.  The  constant  density  case 
shows  a  rather  slow  exponential  energy  decay.  This  energy  decay  is  clearly  the  result 
of  viscous  dissipation.  Both  stratified  cases  show  much  faster  energy  decay  initially, 


Figure  32:  Vertical  energy  flux 


Figure  33:  Total  energy 
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followed  by  an  oscillatory  but  gradual  decay.  This  enhanced  decay  rate  is  caused  by 
i  waves  rapidly  transporting  energy  away  from  the  vortices,  finally  reaching  the  upper 
,  boundary  where  the  energy  ’escapes’  by  being  absorbed  by  the  damping  layer. . 


6  Conclusions 

Several  conclusions  can  be  drawn:  .  'i  • 

1.  A  vortex  pair  in  a  strongly  stratified  fluid  rapidly  loses  energy  to  internal  wave 

motion.  :  .  " 

2.  The  original  vortices  move  horizontally,  instead  of  the  vertical  motion  in  a 

constant  density  flow.  !;i  i: 

3.  The  narrow  region  above  the  interface  between  layers  of  constant  Brunt- Vaisala 
...  frequency  experiences  an  enhanced  level  of  vorticity  during  the  motion.  .  T  ; 
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